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Abstract 

We recall four open problems concerning constructing high-order matrix- 
exponential approximations for the infimum of a spectrally negative Levy 
process (with applications to first-passage/ruin probabilities, the wait- 
ing time distribution in the M/G/1 queue, pricing of barrier options, etc). 

On the way, we provide a new approximation, for the perturbed Cramer- 
Lundberg model, and recall a remarkable family of (not minimal order) ap- 
proximations of Johnson and Taaffe [JT89], which fit an cirbitrcirily high 
number of moments, greatly generalizing the currently used approxima- 
tions of Renyi, De Vylder and Whitt-Ramsay. Obtaining such approxima- 
tions which fit the Laplace transform at infinity as well would be quite useful. 



keywords: Levy process; first passage problem; PoUaczek-Khinchine for- 
mula; method of moments; matrix-exponential function; admissible Pade 
approximation; Johnson- Taaffe approximations; two-point Pade approxima- 
tions 

1. Introduction 

Motivation: The problem of approximating distributions based on em- 
pirical data like moments is one of the bread and butter problems of applied 
probability. 
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In risk theory and related first passage applications (queueing, math- 
ematical finance, ...) there is special interest in approximating densities of 
nonnegative random variables by affine combinations of exponentials, 
also called GHE densities (generalized hyper-exponential), or, more 
generally, by matrix exponential distributions (which allow for "coUision 
of exponents" ) . One reason for that is that this class captures the asymptotic 
behavior in the important "light tails case". 

There exists already a quite extensive literature, based on inverting the 
explicit PoUaczek-Khinchine formula for the Laplace transform, which as- 
sumes however complete knowledge of the model, for example knowledge of 
all its moments. 

Since input data are never certain, it is interesting to develop approxi- 
mations based on finitely many moments (the coefficients of the power 
series expansion of the Laplace transform around 0). 

Some well-known such approximations in risk theory and queueing are the 
Renyi, De Vylder, Gamma, and Whitt -Ramsay approximations, ob- 
tained by fitting one, two or three moments, and the minimal (but arbi- 
trarily high) order three moments fitting formulas of Bobbio, Horvath 
and Telek [BHT05]. Also useful is the Cramer-Lundberg approximation 

*(x) ~ Ce"^^ (1) 

where ^{x) denotes the ruin probabihty (5) and —7 is the so called adjust- 
ment coefficient (i.e. the largest negative root of the Cramer-Lundberg equa- 
tion (13)), which captures the asymptotic behavior in the case of light tail 
claims. All these approximations may be derived from the explicit Laplace 
transform provided by the PoUaczek Khinchine formula (6). 

Producing higher order approximations in the " intermediate regime" when 
a finite, but larger number of moments is known, seems a very challenging 
problem. 

Question 1. In view of the scarcity of approximations fitting more than three 
moments, it is natural to ask what are the difficulties blocking the 
development of high order moments based approximations ? 

Below, we examine this question in the context of first passage theory for 
spectrally negative Levy processes. 

Pade, two-point Pade and other rational approximations. It turns 
out that most of the approximations currently used are Pade approximations 



2 



of Laplace transforms, and that higher order approximations are quite easy to 
obtain (for example using the Mathematica command PadeApproximant), 
since the conversion from moments to a rational Laplace approximation re- 
quires only solving a linear system. It is also quite easy to produce Pade 
approximations with specified limiting behavior of the Laplace transform at 
oo, so called two-point Pade approximations - see Example 2. 

Other rational approximations of interest are those minimizing in least 
squares sense the sum of the coefficients - see Beylkin and Monzon [BM05]. 
Another interesting class are "Tijms approximations", which try to incor- 
porate moments fitting with including the exact Cramer-Lundberg asymp- 
totics (1) as dominant term -see [Wil98]. These may also be obtained by 
using a Pade approximation with a prescribed pole. 

The admissibility of Pade approximations in probability. While 
in principle a great tool due to their easiness of implementation, and their 
convergence for large n, Pade approximations (and variations, like two point 
Pade) applied to Laplace transforms in probability have the drawback of the 
difficulty to guarantee "admissible inverses", i.e. nonnegative densities 
and non-increasing survival functions, when fitting three moments or 
more is desired (note that for fitting two moments mi,m2 of a nonnegative 
random variable, the admissible Gamma approximation 

fix) ~ fix) = '\^/ . i^e a = n = 2 

i ya) 1712 — TTii m2 — m| 

provides an easy solution). 

Even ensuring the nonnegativity of combinations of exponentials 
with fixed given rates is quite a difficult question (since this involves an 
infinite number of constraints), still open nowadays, except for two exponen- 
tials (when nonnegativity of /(O) and nonnegativity of the coefficient of the 
asymptotically dominant exponent are clearly necessary and sufficient), and 
for three exponentials [DL82]. 

Example 1. Consider the example due to Harris [HWR92] (see also [FacOS, 
Ch. 54]) 

fit) = 2e-* - 6e-2* + Ge'^* - 26"* - 3(26"^*) + 2(3e"^*) 

^ Fit) = 2e-* - 3e-^* + 26^=^* 

with canonical coordinates (2, —3,2), and Laplace transform 
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2(5^ + 25 + 3) _ 2 6 1 6 2 3 

~ (s + l)(s + 2)(s + 3) ~ 3(s + l)(s + 2)(s + 3) 3(s + 2)(s + 3) ^ 3s + 3 

S'mce t/ie "Coxian coordinates" (2/3,-1/3,2/3), produced by the partial 
fractions decomposition above are not nonnegative, this is not a phase-type 
distribution of order 3. 

However, by an admissibility criteria for combinations of three nega- 
tive exponentials due to [DL82] ^ (see also [FacOS, BFT08] for a criterion 
which allows colliding exponents) , we know this is a proper density. In fact, 
it is a phase-type density of order 4 - see (2) 

Admissibility by phase- type representations. One "lucky case" in 
which admissibility is automatic is when one has obtained somehow any 
phase-type representation PH{a, A) with a a probability vector and 
A a subgenerator matrix (satisfying Aij > for i ^ j and Al < 0). In that 
case, nonnegativity follows from the probabilistic interpretation of f{t) as the 
density of the absorbtion time of the corresponding Markovian semigroup. 

However, determining when a phase-type representation exists is again 
notoriously difficult, the so called positive realization problem of systems 
theory ^ . 

In example 1, it is possible to show that a phase-type of minimal order 4 is 
available, by a recursive approach of decomposing the admissibility domain 
as union of higher order admissibility polytopes associated to phase-type 



^with exponents equal to i = 1, 2, 3, this is —002 < 2^a;ia;3, or, after normalization 

LOl + 1 < LO3 < UJl + 1 + 2y^UJl. 

^Finding the minimum possible order of such a representation is even harder, and 
known currently only for three moments fitting representations). Attesting further to 
the difficulty of providing admissible phase-type approximations are several interesting 
recent approaches, like the recursive minimal order three moments fitting formulas 
of [BHT05], and the EM algorithm approach [WLll]. 
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representations - see [0'C93].One phase- type representation is: 

a = (1/2, 0,0, 1/2) (2) 



/-II \ 
0-220 
0-33 

V -4 y 



BUTools http:/ /webspn. hit.bme.hu/ telek/tools/butools/butools.html only 
obtains a representation of order 5, further illustrating the difficulty of this 
problem. 

Question 2. Currently, no algorithmic approach for testing admissibility of 
combinations of more than four given negative exponentials (or four terms 
Miintz polynomials) is available ^ . 

The question of providing admissible approximations (with non-specified 
rates) fitting more than four given values of a Laplace transform (for example 
moments) is a priori even more difficult. However, two remarkable excep- 
tions in which this the challenging admissibility problem was solved are the 
minimal order three moments fitting admissible approximations of 
Bobbio, Horvath and Telek [BHT05], obtained by a recursive approach on the 
order, and the non- minimal ones of Johnson and Taaffe [JT89] , an outcome 
of the classical moments theory, which work for any number of moments. 

Contents and contributions. This problem was motivated by the 
desire to provide new admissible matrix exponential approximations in ruin 

theory- 
Necessary ruin theory background, including the Pollaczek Khinchine 
transform, is reviewed in Section 2. Section 3 presents the "key charac- 
ters" in our ruin application, the aggregate loss and its moments - see (16), 
(15). 

First order Fade approximations of the Follaczek Khinchine transform 
due to Renyi and DeVylder are reviewed in Section 4. We also provide 
here a new approximation, for the perturbed Cramer-Lundberg model - see 



§a solution might however be possible by the approach of Faybusovich [Fay02, Thm 4,5] , 
who offers a general representation of the Koszul-Vinberg characteristic function [Seel2] of 
the positivity cone generated by any Chebyshev system, as PfafHan of a matrix of multiple 
integrals (and taking logarithm yields a "self-concordant barrier" function). 
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Theorem 1. The remarkable Johnson-Taaffe approximations are reviewed in 
Section 5, Theorem 2. Using their approach, we may "update" a second 
order approximation due to Ramsay to make it work for arbitrary claims 
having three moments - see Theorem 3. 

In the case of random sums, it is possible to apply the Johnson-Taaffe 
approach both to the individual summands - we call this a " Ramsay- type ap- 
proximation" , and directly to the sum - yielding "Beekman-Bowers-type ap- 
proximations". Theorem 4 in Section 6 provides a comparison between these 
two methods, by comparing their explicit 3 moments Johnson-Taaffe or- 
ders. 

Section 7 discusses two-point Fade approximations -see Theorem 5, whose 
admissibility is an open problem for the moment. Finally, in Section 8 we 
provide some numerical examples. 

2. Ruin theory background 

The perturbed Cramer Lundberg risk process models the reserves 
of an insurance companyby: 



used in collective risk theory to describe the surplus X — {X(t)^t > 0} of 
an insurance company. Here, 

1. uis the initial capital, 

2. ct represents the premium income up to time t, 

3. Zk are i.i.d. positive random variables representing the claims made, 
with cumulative distribution function and density denoted by F[x) and 
/(x), and (some) moments denoted by rrii^i = 1,2, 

4. = {Nx{t),t > 0} is an independent Foisson process with intensity A 
modeling the times at which the claims occur, and 

5. W{t) is an infinite variation spectrally negative perturbation, for ex- 
ample a standard Wiener motion, and o" > is a scale parameter. 

Since the jumps of X are all negative, the moment generating function 
]E|^g«^(*)] exists for all s > and t > 0, and is log-linear in t. The sym- 
bol/Laplace exponent/cumulant generating function k{s) is defined by 





(3) 



K 



(s) = log (Eoe^^(^)) = s{c- XF*{s) + k{s)) , 



(4) 
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where F*{s) denotes the Laplace transform of the survival function of the 
claims, and sk{s) is the symbol of the perturbation aW{t). For example, in 
the Wiener case k{s) = ^s. 

Let T be the first passage time of a stochastic process X{t) below 0: 

T := mi{t > : X(t) < 0}. 

The objects of interest in classical ruin theory are the "finite-time" and "even- 
tual" ruin probabilities 



^{t,u) = Pu[T < t], ^{u) = Pu[T < oo] = Po[L > u 



where L = —2L, also called maximal aggregate loss, is the negative of all-time 
infimum of the process (3), started from 0. 

The ultimate ruin probability "^{u) for (3) is not identically 1 iff the Levy 
drift/profit rate 



p :— c — Ami > := Xmi9 (5) 



is positive, in which case adding the condition lim^^^ \E'(a;) = determines 
it uniquely. 

The Pollaczek-Khinchine formula. Taking Laplace transform of the 
Kolmogorov equation for the ultimate ruin probabilities of the perturbed 
Cramer-Lundberg model yields the Pollaczek-Khinchine formula: 

^' s k(s) s(1- />/;(«) +k(s)/c) ™ 

The first expression, in terms of the symbol k{s) of the Levy process involved, 
is valid for all spectrally negative Levy processes, - see for example [Kyp06]). 
The second emphasizes the fact that the result in the Cramer-Lundberg case 
depends only on the "equilibrium density" of the claims, defined by 

Ux) := F{x)/rm, (7) 

the estimation of which may be a convenient starting point. The third ex- 
pression is equivalent to ^(,x) = ■ip{u)du, x > where ip{u) is the 
density of the continuous part of the distribution of the aggregate 
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loss L (while '^{x) is the survival function of L). The PoUaczek Khinchine 
formula for 



1-p 



(8) 



In the case a = 0, a beautiful probabilistic interpretation of the (non- 
perturbed) Pollaczek Khinchine formula (8) was discovered (independently) 
by Benes, Kendall and Dubordieu, by expanding the denominator into a 
geometric series: 



i-p 



1 - Pfe{s) 



k=0 



(9) 



where fe{x) = F{x)/mi is the equilibrium distribution of Zi. This reveals 
that 'ip*{s) is the Laplace transform of a geometric sum of convolutions of 
the equilibrium/stationary distribution, i.e. ijj{x) = (1 — p) Yl'^=o fe*'''\^) 
^ , which may be visualized by examining the "ladders" of the paths (the 
amounts by which the process —ct + S(t) jumps to new maxima). 

This is the so called Pollaczek Khinchine "ladder decomposition" L = 
X]^=i of the "maximal aggregate loss random variable" L, with N a ge- 



ometric r.v. Pr[A^ = k] 
number of ladders. 



:i-p)p\k = o,i,...,p 



Ami 
c 



, representing the 



Remark 1. The ladder decomposition, cr > 0. In the Brownian per- 
turbed case, a beautiful probabilistic interpretation of the Pollaczek Khinchine 
formula (8) was recently discovered by Dufresne-Gerber and rederived in an 
elementary way by Kella, by rewriting (8) as: 



1 



i-p 



1 - pftis) + h{s)/c 1 + k{s)/c 1 - p/*(s)t^ 



(10) 



s)/c 



reflects the fact that L is an independent sum of a "first creep at the current 
infimum" (which in the diffusion case is an exponential of rate |§j and of an 



§ Another interpretation of i/'*(s) is as Laplace transform of the stationary waiting time 
of the M/G/1 queue -see for example [Asm03, Thm VIII. 5. 7]. 
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alternating geometric sum of "compound Poisson ladders and further creeps" 
when a>0 - see [DG91, Fig 2], [KelllJ. 

In the case that the Brownian perturbation aW{t) is replaced by a general 
spectrally negative perturbation Y with non-zero expectation E[Y{1)] > 
0, a similar ladder decomposition holds true. Let 

K*{s) — Sk{s) 

denote the Laplace exponent of Y. Note that by the Wiener-Hopf factorisation- 
see e.g. [Kyp06], k{s) is Laplace exponent of the possibly killed downward lad- 
der process (and s is the Laplace exponent of the up-crossing ladder process) , 
that is, k takes the form 

^s)^c + Ss+ {l-e~^'')u{dx), 

J{0,oo) 

where the killing rate c and the drift 5 are non negative constants and the 
Levy measure v satisfies the integrability condition 



/ [1 A x]i'{dx) < oo. 

J(O.oo) 



Then, the formula (10) still holds, providing a decomposition of L as an 
independent sum of the increment of the ladder process of the pertur- 
bation Y at an independent exp(c)- random time and the geometric 
sum of further such increments and "compound Poisson ladder height incre- 
ments". In particular, if the perturbation is given by a completely asymmetric 
stable process, i.e. k*(s) = s" , a G (1,2), then k*_{s) = s°^^ and we identify 
i+K*^(s)/c '^^ Laplace transform of the non-negative random variable ^^(c)^' 
where Y°'~^ is a stable subordinator with parameter a — 1. 

Remark 2. Note that the aggregate loss L is the mixiure of a discrete mass 

of i/j* (oo) = 1 — ^1/(0) = 1 — p at 0, and of a continuous random variable. 
Letting i(j*{s) denote the Laplace transform of the density 4>{u), note the 
decomposition 

i/j^s) = 1 - p + pip*{s) <^ i/jix) = (1 - p)5o{x) + pi/jix) 

r(s) ■= = ns)-{i-p) f:{s) 

^*(o)-^*(oo) p ^h-pft{sy^''^ 
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where we denoted by ip{u) the inverse Laplace transform ofip*{s), given by 



The behavior at oo of ip*{s) distinguishes between the nonperturbed (a = 
0) and perturbed case (a > 0): 



The Laplace transform ip*{s) is an essential quantity in the theory of Levy 
processes 

Remark 3. The roots of the denominator k,{s) = in the Pollaczek Khin- 
chine formula determine the asymptotic behavior of ultimate ruin probabili- 
ties. More generally, an important role is played by the roots of the Cramer 
Lundberg equation 



Remark 4. In practice, the true distribution of the claims ( and interarrival 
times) is of course unknown, and since the Pollaczek Khinchine formula 

requires this knowledge, it should be viewed more as a theoretical than a prac- 
tical tool. It may be argued that the most reliable information available in 
insurance data is contained in the first few integer moments, and thus it seems 
natural to start building approximations by fitting moments, or, equivalently, 
by Fade interpolation of the Laplace transform at the origin. 

3. The moments of the aggregate loss 

Prom now on, we will assume a classic Brownian perturbation. As an 

alternative to Ramsay's approximation of the ladder distribution, wc may 
approximate directly the aggregate loss distribution. One advantage is that 
the second moment of the aggregate loss satisfies automatically the second 
order represent ability constraint Cy > 1/2. This allows focusing on the third 
moment constraint. 





(12) 



k{s) — q,q > 0. 



(13) 
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Consider the expansion 

^*{s)= + - ^ (14) 

r{s) ^ 



P + V2,aS/2 - 773SV6 + 

e 



V ^mi 2\mi' I mi ' mi ' / 
^ W 

1 - p (1 - + 2&)V2 + ^«V6 - ::^«V4! + ...) 

where rji — Xrui^i — 0,1,2,... are the moments of the Levy measure, and 
?72,a = 772 + o-^- 

Remark 5. When (7 = 0, the expression in the last parenthesis of (14) 

/las moments Thi — , identifying Li as the famous equilibrium/ stationary 

excess/ladder variable generated by Zi, with density fe{x) — F{x)/mi, and 
stationary excess moments fhi = ^ . 

Remark 6. Let Li = EL*, Aj = li/i\,i > 1 denote the moments and "facto- 
rially reduced moments" obtained by "normalizing" with respect to the expo- 
nential moments of the maximal aggregate loss. These may be easily obtained, 
either by the recursion equivalent of the Pollaczek Khinchine formula: 

n-l 

a\ ~ I ~ \ ^ 1 ~ "^'^ 

UAn — //n + / ^ IJ'k^n-k, n- ^ L, IJ,k ■— -jy — 

fc=l ^ 

or by expanding (14) in power series. The first factorially reduced moments 
are: 

A.=EL=!|l + ^ = ^!^, (15) 
Bp 2p 

X,-EL /2!-^ + \, A3 -EL /3!-^ + 2^^ + ^ 

Tn,r4/., "^4 , 2mim3/3! + (m2/2)^ 3m?7ft2/2 mf 
A4 = EL /4! = ^ + + ^3 + ^' 

an(i t/ie mass 0/ t/ie continuous part zs Aq := Pr[L > 0] = p = 
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Remark 7. When a — 0, the factorially reduced moments admit an in- 
teresting decomposition: 

^ oo 

ris) = (1 - ph — = 5^(1 - p)pv:is)r = 

1 - Pfe{s) k=0 
oo oo / 

k=l k=l ii+i2+---=k 

where Y^' denotes sum over all decompositions of k as a sum, and 

+ 1)! z! 
where m := (^^^, := ^ anc? /ii := ^. 

Remark 8. The moments of the conditioned continuous r.v. L\L > 0, nec- 
essary for applying certain results from the literature, may be obtained by 
dividing by p — . 

Remark 9. Note that the corresponding moments li may also be viewed as 
moments of the ruin function: 

poo poo ~ ~ 

y^fe :- Afc+1, ... 

4. Renyi, De Vylder, and a new simple approximation for ruin 
probabilities 

Approximations of ultimate ruin probabilities. The problem of 

approximating ultimate ruin probabilities ^(m) for the Cramer Lundberg 
model (3) using data on the distribution F{u) of the claims is a classic of 
applied probability, dating back before 1900. Its roots may be traced back 
to the Danish mathematician TN. Thiele, who founded the first insurance 
company, Hafnia, who is also the inventor of cumulants and of Thiele 
continued fractions useful for Laplace transform inversion [Coh07] . 
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In this section, we make the observation that the exponential mixture 
approximations recalled in the introduction are particular cases of Pade ap- 
proximations of Laplace transforms (a theme already present in Thiele's re- 
search preoccupations). We also provide a new simple approximation in this 
vein for the perturbed Cramer- Lundberg model in Theorem 1. 

Remark 10. Concerning our application which involves random sums, we 
make the observation that it is possible to use Pade, Johnson-Taaffe and any 
other rational approximations of Laplace transforms at three levels: 

1. for the density of the claims, based on the estimates ofmi,i > 1 

2. for the equilibrium density of the claims, using the estimates of the 
equilibrium moments Thi,i > 1, and the profit rate p. The second level 
is intuitively superior to the first, since the equilibrium density is mono- 
tonically decreasing, even when the claims density isn't ^ . We will call 
this Ramsay type approximation. 

3. for the aggregate loss density transform ip*{s) (or, equivalently, the ccdf 
transform \E'*(s)j, using directly the moments Aj+i,i > 0, given in (15). 
This amounts to working directly with the Pollaczek-Khinchine 
formula, instead of approximating its denominator, and so intuitively, 
should be better, at least under certain conditions. ^ . We will call this 
Beekman-Bowers type approximation. 

One Pade approximation we consider here is: 



where Vade{jn,n) denotes the "classic Pade approximation" based on the Tay- 
lor series around 0, where m = n or m = n -|- 1 for the classic/perturbed 
Cramer Lundberg process, respectively, and where [/(s)]Ar denotes trunca- 
tion of a power series to its first terms, with N = m -\- n -\- 1 (for "theo- 
retical models" where an expression for 'il)*{s) is available, we may also take 



^For example, this gives rise to a smaller JT index, as illustrated in Example 5 
^More precisely, the third method is expected to be better in the case of light tails 
claims and heavy traffic, while the second is expected to be better in the case of heavy tails 
claims and light traffic. The fact that both methods are better sometimes is illustrated in 
Theorem 4, from the point of view of yielding a smaller JT index. Identifying the domains 
within which methods two and three are preferable in "boundary cases" is not an easy 
task. Note also that a mixture of the two has been also proposed [Sak04]. 




p 



(16) 
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N — oo). The first case is applicable to the classic, and the second to the 
perturbed Cramer Lundbcrg model. As mentioned, the motivation of (16) 
is that the classic De Vylder approximation, is precisely the one point Pade 
approximation of around 0, of orders (n — l,n), with n = 1. 

A second class of Pade approximations we experiment with is: 

{lp*)^''\s) = 'Pade^^rn,n,mi,m)){'^*{s)) = 'Pade^(^rn,n,mi,m)){[ / N / ]iv) (17) 

K[S J / S 

where 'Pade(^(^m,n,mi,ni)) denotes a two point Pade approximation based on the 
power series around and oo. These are indispensable when dealing with 
the perturbed model. 

Example 2. The simplest Ramsay-type approximation is the one moment 
Renyi exponential approximation of the equilibrium density (which 
may also be viewed as a Pade (0, 1) approximation of the aggregate loss 
density, which imposes also the correct limiting behavior limg^^ s'if* (s) — 
^(0) = p — 1/(1 + ^) of the Laplace transform at oo). This amounts to 
looking for an approximation of the form 



s + bo ' s + 6o/(l -p)' 

where f*{s) denotes the Laplace transform of the stationary excess density 
of the claims (7) (the two being related by the Pollaczek Khinchine formula 

Fitting one moment yields ^ = mi, where rhi — (the first moment 
of the stationary excess density), and 

**(s) ~ 13/ — , ^ -^{x)^ pe-''^^-P^'^\ (18) 
~ ~\~ s 

mi 

which is also correct in the limit p — >■ 1 when rhi exists see [Kal97], [GraOO, 
(31)]. 



De Vylder's exponential approximation. The simplest Beekman- 
Bowers type approximation is De Vylder's (19), one of the most popular 
approximations for ruin probabilities, due to its simplicity and asymptotic 
correctness [GraOO], despite of its its being expressed in terms of only the 
first three moments of the claims. 
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Derivation of De Vylder as a two moments Pade approxima- 
tion. Dc Vyldcr's approximation was obtained originally by equating the 
first three moments of the original Cramer Lundberg process with those of 
a new process with exponential claims, and with different arrival intensity 
and premium rate. We check now, using the moments (15) of the aggregate 
loss density, that De Vylder 's formula coincides with the Pade (0, 1) approx- 
imation around of the Laplace transform ^*(s) of the ruin probabilities. 
We start by expanding in power series the numerator and denominator of 
the PoUaczek-Khinchine formula (6): 



Am2/2-Am3s/6 + ... _ a 



p + Xm2s/2 — + ... s + a 

<^ as{p + Am2s/2 — Amas^/G + ...) 
(s + a){Xm2s/2 — Xm^s'^/Q + ...) 

ap = Q;Am2/2, 
am2/2 = 1712/2 — am^/d 




3Am| 
3Am2+2pm3 ' 
6pm2 



3Ams+2pm3 



and so 



3Am^ 



'''' ^ + « «+3A4f|^ s(3Xml + 2pms) + 6pm2 ^ ' 

Question 3. Identify domains within which methods two and three are "prefer- 
able" in some sense, in the simplest case of exponential approximations, i.e. 
compare the Renyi and De Vylder approximations. 

Example 3. The Cramer Lundberg process with Brownian pertur- 
bation aB{t). The Laplace exponent is 

<s)^ps + ^s^ + J2i-sr'^, (20) 

k=3 

where Xruk are the moments of the Levy measure, and K2 — Xm2 + cr^. 

Now we have two further unknowns of interest, the probability of "creeping 
ruin" and that of "ruin by jump". The respective Laplace transforms satisfy 

^^(^) = (^) = (^)' = ^ - ^ - (^) (21) 
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The Pade K = 0, L = 1 approximations are unreasonable, since they 
cannot satisfy the boundary conditions ^d(O) = 1, ^j(O) = 0. To satisfy those 
as well as the equation 

^}{s) = "^ris) ^ 1 - s^*{s) = {s+ f^ms) (22) 

which follows from (21), we must use at least a K — 1, L — 2 approximation. 
Theorem 1. Consider the exponential approximation 

^a{x) = ^l^e-.- + fl^e-^^ (23) 

S2 — Si S2 — Si 

m.[x) = % ]^^-S,X _ ^-S2X^^ (24) 

S2 — Si 

for the "creeping ruin" and "ruin by jump", which satisfy ^'^(0) = 1, ^'j(O) = 
0. Then, by fitting the first two moments of the aggregate loss L, one is led 
to the following admissible approximation: 

3m2 Am2 3Am| 

and —si and —S2 the negative roots of + (a^ + aj + ||)s + ad|§ = (whose 
discriminant is non-negative). 

Proof: \l/*(s), ^rf(s) = s2^^^f_^i,Q ^'^^ quotients of manic polynomials (to 
satisfy liiiis^^oo ^^^(s) = ^(i(0) = I), with three free coefficients, but the 
second condition in (22) imposes one more condition bo = (so that 

"if^s) — ^ — j^), leaving only two free coefficients. Finally, fitting the 

cr 

first two coefficients around of 

2p X 2p s + 

^ '^'^''^s2 + (a, + a, + J)s + a,5 

yields: 



2p \ ad' 2p ad o'^ 

(T^ 0-j Ama 

2p 6p 

with solution (25). Then, assuming si < S2, Laplace inversion yields (23). 
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Moreover, both and are admissible. Indeed, this is obvious for "^j, 
since its initial value and its dominant coefficient are non-neqative. 

The same is true for \E'rf; indeed, we may check that its dominant co- 
efficient is non-negative, i.e. that si < a^, by noting that s^ + (a^ + aj + 
2p/(j^)s+a<i2p/(j^ is negative ats — —a^ (since (ad)^+(ad+aj+2p/(T^)(— 0^) + 
a^'pja^ — —adUj < 0, with and Uj being positive). Therefore, Si < < 

S2- 

Remark 11. It is easy to check that this approximation is exact for exponen- 
tial claims. Indeed, in that case the density transform in the parametrization 
(26) is: 

s + m^^ 
(s + mr')(l + gs) + ^smi 

2p s + m]"^ 

" + s{^m, + fm-,' + 1) + mr^ ^ + + + |§) + mr^|§ 

with Qfi — l/mi,aj — 2Xmi/a'^, and it is easy to check that this coincides 
with our approximation (25), in the case of exponential claims. 

Remark 12. Another admissible approximation exact in the exponential 
case, but fitting now only one moment is aj — 2Xmi/a^, aa — 1/mi = 
(with similar admissibility proof). 



P 



1 + ^s + ^s^^ 
S + m'^^ 



5. Johnson- TaetfFe approximations 



This elegant moment fitting method approximates by common order 
Erlang mixtures approximations ^ : 



^ poo 

fi^) ~ w;ie„(pi,x) := / en{p,x)uK 
i=i -^0 



(dp) 



(27) 



^One hint which may explain this choice comes from discrctizing the well-behaved 
Post-Widder Laplace inversion method [ Jag78] , which leads to Erlang mixtures of common 
order. Note also that pure Erlang distributions are the unique extremal points minimizing 
the coefficient of variability, within phase-type distributions representable at order n, and 
that this is a dense subclass of phase-type distributions. 
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where en{p,x) is the Erlang density. Finding a discrete measure vxidp) with 
positive coefficients Wi requires solving a classic De Prony system with respect 
to Erlang reduced moments 

/^i = — , = — — TT, Pk = , n^k) = n{n + l)...{n + k - 1) 

n n[n + ij '/^(fc) 

obtained dividing by the moments n(^k) of the scaled Erlang density of shape 
n. 

Example 4. When n — 1, one finds the factorially reduced moments 

mi mi mi 

which intervene in constructing GHE approximations (which have positive 
weights, if the moments m^ are far enough from the boundaries of the Stieltjes 
space of moments > 1, >l). 

It turns out that Erlang reducing leads to "positive moment sequences" 
Pk satisfying ^ > 1,^^^ > 1,.... (i.e. having positive Hankel matrices of 

moments up to any desired degree), for n big enough. Then, the classic 
Stieltjes method of creating discrete moment fitting measures with positive 
weights may be apphed. 

Theorem 2. Johnson- Taaffe three moments approximation. Letm — 
(mi,m2,m3) denote the first three moments of a nonnegative r.v. Let 

n* = fmax{(mom2/mi — {mim^/ml — — 1}] 

1 2m2 - ^3 -I 

= maxj^^ —}\ 

m2 - 1 m3 - m2 

denote the explicit " JT index of degree" 3 giving the smallest "JT sector" 
- see [BHT05, Fig 3-7] - containing our target moments. Let /j, — {pi,i — 
1, ...3) denote the Erlang reduced moments of order n* , and let 



bo = P2- Pi, bi= p3- P1P2, &2 = P3Pi - P2 



denote the coefficients of the denominator of the classic second order Fade 
approximant, and let 

2bo 1 - sipi 

Sl,2 = ,^ ^, , , = :<^2 = 1-^1 

bi If ^ybl - 46o&2 *2 - Si 
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denote its roots and ( non-negative ) partial fractions coefficients. 

Then, the Erlang reduced moments satisfy /j, e {60 > 0, 62 > 0}, and 

2 



is a nonnegative, decreasing "common order Erlang n* " density fitting our 
three moments. 

Using this result, we can "fix" Ramsay's approximation by sums of two 
exponentials [Ram92] to make it work for any valid three moments, at the 
price of using higher (not minimal) order approximations - see Theorem 3. 

Remark 13. In principle, 2k + 1 moments fitting Johnson-Taaffe 
Pade approximations could be obtained by determining numerically higher 
order JT indices ensuring the positivity of corresponding Hankel determi- 
nants. The Hankel determinants that must be taken into consideration are: 

■MkiV) = {7724(7722 — ml) + 7712(77737711 — 7712) + 7773(77127711 — 7713) > 0, 
m^nii (77117713 — 77I2) + {—ml + 2m2m4m^ — 771177I4) > 0, } (28) 

The inequalities above define implicitly the JT index of degree 5. 

Theorem 3. Ramsay updated Consider the ruin problem for the classic 
(nonperturbed) Cramer- Lundberg model (3). 

1. Imposing the correct limiting behavior liuis-^^ s'i>*{s) — ^(0) — p — 
1/(1 + ^) at 00, the Pade (1,2) approximations for the equilibrium 




i=l 



density and ruin transforms are: 




N 



bo + ais 



D 




D-N 



62S + 61 - ai 



P 



D-pN 



P 



625^ + biS + bo' 



bo = p2-Pi, bi = p3-p2Pi, 
ai = 61 - pibo = P3- 2/ii/i2 + pI 
bis = 61 - pai, 60 = (1 - p)bo, 



b2 = P1P3 - P2 



where p2 are factorially reduced equilibrium moments. 
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2. For the approximate inverse transform fe{x) to be hyperexponential (in 
particular, nonnegative and completely monotone), it is sufficient that 
bo > 0,62 > 0. 

3. Higher order moments fitting is possible by Erlang approximations of 



ordern* biq enough, where for three moments n* = ["maxj '"^ 2 , — — — j- 

1)1 • 



6. Comparison between the JT indices of Ramsay and Beekman- 
Bowers type approximations 

In this section we present a simple comparison between these two ap- 
proaches, based on comparing their Johnson-Taaffe indices. 
The moments of the aggregate loss L satisfy 



r2 = !? = 2+^^^ 



3+ - 



2 + > 2 + ^ > 2, 



21711 _i_ Q/mi \J 

e2 + e J 



3 + 



2 + ^7712 



(29) 



The second moment satisfies the necessary inequality ^2 > ^ [BHT05] 
already with n = 1. We only need to investigate the necessary inequality for 
the third moment, which is: 



n 



> 



n+2 
n+1 



J(L) := 



2/2 - I3 

I -To 



— 9 



^2 - im,-,0 - 2 



^7712 — ?Tl2m3^ — 2 



< 




9'^ \mo — m2m3 
01712 ( Orhs + 4^+6 



[e^2 + 2) 



Here, J(L) — is a "partial JT index", based on the third moment 

I3—I2 ^ 

admissibility condition, and the "normalized moment" 77(L) = ^ is a mono- 
tone transformation (since J = ^5? is a decreasing function in the relevant 
range 77 e (1, 00)), which has been already used in the literature. 

A rough indication of the performance of the Beekman-Bowers and Ram- 
say methods will be obtained now by checking which of r]{L), r){Li) is higher. 
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Theorem 4. a) The partial J index of the aggregate loss is strictly smaller 
than J(Li) = iff 

01712 < 



3/2^2 - ^3 3/2 - r]{Li) 



— m2 



v{Li) - 1 



In particular, if ri{Li) = ^ > ^, in which case ri{L) < | holds as well, a 
two terms exponential mixture distribution matching the first three moments 
exists only for the equilibrium ladders Lj. 

b) More generally, letn denote the unique integer such that < r]{Li) < 

> 2, is satisfied. If furthermore Xi < 6m2 < ^2, where Xi^2 — 



^?H-i)(>-(n+2) — ^' then a n'th order Erlang mixture distribution matching 
the first three moments exists only for the equilibrium ladders Lj. 

c) For any 9, there exists a mixture of Erlang n distribution matching the 
first three moments of L if 

n > —. 



mz — m2 

Proof: a) Consider the unimodal function 

U ax'^ + 4x + 6 



r}{a,x) ^ ^ ^ 



2 



l2 (^ + 2) 

where a = ^,x = 6*7712, which takes values 
3fl 2 

7: T < ^(o, < max[77(a, 00) = a, r]{a, 0) = 3/2] 

(the lower minimal value is achieved for x* = -x^ 6* — ^ ^ ^ ). 

b) Fitting L is possible when 

ax^ + 4x + 6 n + 2 ,, . , n , „, , 

^ > ^ ((^ + 1)« - + 2))^ - 4^ + - 1) > 0. 

In the prescribed range of a, the discriminant A = 2(n^ + n — a{n'^ — 1)) 
is always positive, and the inequality holds when x < xi when x > X2, 

, 2=FV2A/n2+n-a(n2-l) 

where Xi 9 = r^-p; — r—^\ • 

-^'^ (n+l)o— (n+2) 

c) Follows by minimizing r][9) in 9, which yields 9* — ^ and r]{9*) = 

3m3— 2m2 7(0*) = ms 
2m3— 7712 ' ms— m2 
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0.25 
0.20 
0.15 
0.10 
0.05 



0.2 0.4 0.6 0.8 1.0 

Figure 1: Ruin probabilities for uniform claims. Red line: original calculated by numerical 
inverse Laplace; green line: by 3 moments [JT89] on the uniform itself; blue line: by 
3 moments [BHT05] on the uniform itself; yellow line: by 3 moments [JT89] on the 
equilibrium of the uniform; black line: by 3 moments [BHT05] on the equilibrium of 
the uniform. 

Example 5. Consider a three moments fitting example: the U[0, 1] rv, with 
moments mi = \,m2 = ^,"013 = |,m4 = |. The Johnson- Taaffe index is 7, 
and so this procedure yields 14 phases, while the Bobbio, Horvath and Telek 
[BHT05] and He-Zhang methods [HZ06] yield only 9 phases. 

For the equilibrium distribution, the moments are rhi = |,?7i2 = ^,^3 = 

the normalized moments are 

o m2 3 o 9 
m| 2 m2mi 5 

and the Johns on- Taaffe index is 4. The Bobbio, Horvath and Telek method 
yields the minimal order of 3 phases. 

Finally consider the JT approach for a "Beekman- Bowers" approxima- 
tion. Computing the partial JT index, we find 

4 (9^2 + + 10) 
9^2 + 40 ■ 

By Theorem 4 a), this is strictly less than the 4 required for the Ramsay 
approach when 9 <1. Taking into account the integer part, we find out that 
< 1.5 still yields a J index less or equal to 4. 
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7. Admissible two point Pade approximations 

We consider now briefly including further information beyond moments, 
for example by using two-point Pade approximations (of special interest in 
the case of Levy processes with infinite variation paths - see Example ). 
Ensuring admissibility in this case is however an open problem. 

We may formally expand the ruin transform at infinity as well: 



„oo °o k °° 

Jo u-n i.-n 



(30) 



k=0 k=0 

Here, 

*(0) 



1 if (7 > 

ni^p if (J = 

is well known, and the derivatives at 0: 



*'(O) = -^(l-^) = -^ = -^0^, (31) 
c c mi 

*"(0) = -*'(0)(/(0) - ^), 

c 

^(3)(o) = -^'(0) (/'(O) + 2%(0) - (^)2) , 

*(4)(o) = -^'(0) ( f"(0) + 2^/'(0) - ^/(O)' + 3(^)V(0) 
V c c c 



' C ' 



may be obtained recursively, by differentiating the integro-differential equa- 
tion for ^!{x). 

Theorem 5. Two-point Pade-Ramsay approximation Imposing both 
the correct limiting behavior lims^oo s^*(s) = ^I'(O) = p = 1/ {1 + 6) at oo, 
and the first derivative at oo ^'{0) = — ^^^^ leads to the two-point Pade 
(1,2) approximation: 

■''^ ^ bo + bis + b^s-'' ^bo + bis + b2S^ ^ ' 



where o-i = ^ o,nd the coefficients bi,bi, obtained fitting the first two equilib- 
rium/aggregate loss moments, are: 

62 = (2mim3 — 3m2)/6, 61 = {m^ — 3mim2)/3, bg = m2 — 2m^, 
bi = bi- pai, 60 = (1 - p)bo. 
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Proof: The theorem follows from the identity 



s^[s) Pb,s^ + ih^pa,)s + bo{l-p) ^^^^ 

by taking into account that we are looking for an approximation of the form 

**(s) a; p = ^ := p^, 

hs^ + hs + bo D 



which will satisfy lim.^oo = 1, lim.^oo s(sf - 1) = *'(0)/p = -i^. 



8. Numerical results 

8. 1 . Without perturbation 

In the following we illustrate the application of the approximations pro- 
vided in Theorem 3 and 5. 

8.1.1. Mixed exponential claim distribution 
With mixed exponential claim distribution 

_ 315e-5^ 7e-^^ 27e-3^ 3e-2^ Te"^ 
•^^^^ ~ ^28~ ^ ~8~ ^ + + 128" 

and A = l,c = 2/5 the ruin probability is 

^ _ 245e-9^/2 i35e-W2 seTe-^^/^ 735e-3^/2 i9845e-^/2 
~ 32768 ^ 8192 ^ 16384 ^ 8192 ^ 32768 

Wc approximated the ruin probabilities by the Renyi and Dc Vylder's 
first order exponential approximation given in (18) and (19), and by the 
approximations provided in Theorem 3 and 5. In Figure 2 we show the ruin 
probabilities with the different approximations and in Figure 3 the relative 
error with respect to the exact solution. As expected, the approximation in 
Theorem 5 is better near a; = 0, and that may be exploited to obtain a better 
approximation by switching between the formulas in Theorem 5, Theorem 3 
when they cross, starting with the first. 
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Figure 2: Ruin probabilities with mixed exponential claims. 
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0.0001 
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Figure 3; Absolute relative error of ruin probabilities with mixed exponential claims. 
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X 




Rcnyi 


uc V yiQcr 


TlioorciTi 3 


Tlioorcm 5 







QOQOQl 

U. t/Ut/Wt? -L 




QOQOQl 


QOQOQl 


300. 


0.521143 


0.529743 


0.522539 


0.521107 


0.522526 


600. 


0.308668 


0.30869 


0.309273 


0.308713 


0.309268 


900. 


0.182866 


0.179879 


0.183048 


0.182888 


0.183047 


1200. 


0.108338 


0.104818 


0.10834 


0.108347 


0.10834 


1500. 


0.0641841 


0.0610794 


0.0641226 


0.0641869 


0.0641233 


1800. 


0.0380254 


0.035592 


0.037952 


0.0380257 


0.0379527 


2100. 


0.0225279 


0.0207401 


0.0224625 


0.0225272 


0.0224631 


2400. 


0.0133465 


0.0120856 


0.0132948 


0.0133456 


0.0132953 


2700. 


0.00790706 


0.00704247 


0.00786872 


0.0079062 


0.00786908 


3000. 


0.00468448 


0.00410377 


0.00465722 


0.0046838 


0.00465748 



Table 1: Ruin probabilities with Gamma claims (a = 0.01, = 100) and its approxima- 
tions 



8.1.2. Gamma distribution with a = 0.01,(3 = 100 
Let us consider now the Gamma claim distribution 

with parameters a = 0.01, /3 = 100 and with claim arrival intensity A = 
1 and loading factor 9 = 0.1, which appeared frequently in the literature 
([Ram92], [GraOO]). As for the previous example, we approximated the ruin 
probabilities by the formulas given in Renyi, (19) and Theorem 3 and 5. Ruin 
probabilities themselves are given in Table 1 while in Table 2 we provide 
absolute relative errors. 

8.1.3. Gamma distribution with a = 2.5, P = 1 

Next we consider Gamma distributed claims with a = 2.5, /3 = 1 and 
A = 2/5, c = I (-1 + 4^2) . This case is interesting as it has been shown in 
[ACHll] that direct, moment based Fade approximation of the claim distri- 
bution does not result in valid distributions. The approximations presented 
in this paper leads instead to valid ruin probabilities. Ruin probabilities 
themselves are given in Table 3 while in Table 4 we provide absolute relative 
errors. 

8.2. With perturbation 

In this section we illustrate the application of the two approximations 
given in Theorem 1 and Remark 12, respectively. 
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X 


Renyi 


DeVylder 


Theorem 3 


Theorem 5 










0.0288462 










300 


0, 


.0165011 


0.00267814 


0.0000688428 


0, 


.00265325 


600 


0, 


.0000714085 


0.0019599 


0.000146799 


0, 


.00194387 


900 


0, 


.0163373 


0.000993297 


0.000119325 


0, 


.000986105 


1200 


0, 


.0324864 


0.0000177571 


0.0000819863 


0, 


.000019394 


1500 


0, 


.0483709 


0.000957418 


0.0000440626 


0, 


.00094697 


1800 


0, 


.0639946 


0.00193166 


6.12674 X lO"'^ 


0, 


.00191241 


2100 


0, 


.0793618 


0.00290493 


0.000031798 


0, 


.00287691 


2400 


0, 


.0944767 


0.00387725 


0.0000697035 


0, 


.00384047 


2700 


0, 


.109343 


0.00484863 


0.000107635 


0, 


.00480311 


3000 


0, 


.123966 


0.00581907 


0.000145554 


0, 


.00576482 



Tabic 2: Relative error of approximate ruin probabilities with Gamma claims (a = 
0.01, /3 = 100) 



X 


exact ^'(a;) 


Rcnyi 


DeVylder 


Theorem 3 


Theorem 5 


0. 


0.268422 


0.268422 


0.299749 


0.268422 


0.268422 


0.5 


0.22854 


0.217791 


0.237348 


0.22894 


0.228126 


1. 


0.189678 


0.176711 


0.187938 


0.189655 


0.189069 


1.5 


0.154441 


0.143379 


0.148813 


0.154172 


0.154016 


2. 


0.124037 


0.116334 


0.117834 


0.123743 


0.123926 


2.5 


0.0986589 


0.0943911 


0.0933036 


0.0984496 


0.0988216 


3. 


0.0779451 


0.0765868 


0.07388 


0.0778418 


0.0782763 


3.5 


0.0612929 


0.0621407 


0.0584999 


0.0612758 


0.0616894 


4. 


0.0480435 


0.0504196 


0.0463215 


0.0480817 


0.04843 


4.5 


0.0375759 


0.0409093 


0.0366785 


0.0376414 


0.0379079 


5. 


0.0293456 


0.0331929 


0.0290429 


0.0294185 


0.0296037 



Table 3: Ruin probabilities with Gamma claims (a = 2.5, /3 = 1) and its approximations 
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X 


Renyi 


DcVylder 


Xhcorem 3 


Xlicorem 5 




n 










0.5 


0.0470337 


0.0385392 


0.00175079 


0.0018135 


1. 


0.0683674 


0.00917815 


0.000125885 


0.00321361 


1.5 


0.0716262 


0.0364389 


0.00174365 


0.00275474 


2. 


0.062096 


0.0500072 


0.00236637 


0.000891164 


2.5 


0.0432581 


0.0542808 


0.00212105 


0.00164952 


3. 


0.0174272 


0.0521542 


0.00132515 


0.00424903 


3.5 


0.0138333 


0.045568 


0.000278168 


0.00647046 


4. 


0.0494556 


0.0358424 


0.000794076 


0.00804433 


4.5 


0.0887094 


0.023884 


0.00174297 


0.00883326 


5. 


0.1311 


0.0103171 


0.0024838 


0.00879521 



Table 4: Relative error of approximate ruin probabilities with Gamma claims (a = 2.5, f3 = 
1) 



8.2.1. Mixed exponential claim distribution 

As without perturbation, in this case the exact ruin probabihties can be 
calculated by symbolic inversion of the Laplace transform. We applied the 
approximations given in Theorem 1 and Remark 12 for three values of cr, 
namely, 0.1, 0.5 and 2. Figures 4-6 show the exact and approximate ruin 
probabilities and we depicted the two components of the ruin probability 
(by diffusion and by jump) as well. The two approximations result in distin- 
guishable curves only in the case a = 0.1. The relative errors for the three 
values of a are provided in Figure 7-9. The error is smaller for larger values 
of (T and for most values of x the better approximation of ^(x) is by the 
approach of Theorem 1. 

Acknowledgement: We thank Jiandong Ren for useful comments. 
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Figure 4: Ruin probabilities with mixed exponential claim distribution with a = 0.1, with 
the approximation given in Theorem 1 (left) and the one given in Remark 12 (right) 




Figure 5: Ruin probabilities with mixed exponential claim distribution with a = 0.5, with 
the approximations of Theorem 1 and Remark 12 (the two approximations look the same) 
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Figure 6: Ruin probabilities with mixed exponential claim distribution with cr = 2, with 
the approximations of Theorem 1 and Remark 12 (the two approximations look the same) 




Figure 7: Absolute relative error of the approximations given in Theorem 1 and Remark 12 
with mixed exponential claims and cr = 0.1. 
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Figure 8: Absolute relative error of the approximations given in Theorem 1 and Remark 12 
with mixed exponential claims and a = 0.5. 




Figure 9: Absolute relative error of the approximations given in Theorem 1 and Remark 12 
with mixed exponential claims and a = 2. 
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